



IN KIA 



INSTITUT NATIONAL DE RECHERCHE EN INFORMATIQUE ET EN AUTOMATIQUE 



c 

(N 

> 

O 



Formal Proof of a Wave Equation Resolution 
Scheme: the Method Error 

Sylvie Boldo — Fran9ois Clement — Jean-Christophe Filliatre — Micaela Mayero 

Guillaume Melquiond — Pierre Weis 



(N 
> 

(N 
OC 

q 

in 

c 
c 

> 

X 

\- 



N° 7181 

Janvier 2010 




o 

z 

LU 

+ 









05 

to 



CO 



INSTITUT NATIONAL 

DE RECHERCHE 

EN INFORMATIQUE 

ET EN AUTOMATIQUE 




INKIA 



centre de recherche 

SACLAY - IlE-DE-FRANCE 



Formal Proof of a Wave Equation Resolution Scheme: the 

Method Error 

Sylvie BoldqjJ , Francois Clemento , Jean-Christophe Filliatre^* , Micaela 
Mayerdi^l , Guillaume Melquiond*^ , Pierre Weis^ 

Theme : Programmation, verification ct preuvcs 

Observation et modelisation pour les sciences de I'environnement 

Equipes-Projets ProVal et Estime 

Rapport de recherche n° 7181 — Janvier 2010 — [TS] pages 



Abstract: Popular finite difference numerical schemes for the resolution of the one-dimensional 
acoustic wave equation are well-known to be convergent. We present a comprehensive formaliza- 
tion of the simplest scheme and formally prove its convergence in Coq. The main difficulties lie in 
the proper definition of asymptotic behaviors and the implicit way they are handled in the math- 
ematical pen-and-paper proofs. To our knowledge, this is the first time this kind of mathematical 
proof is machine- checked. 
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Preuve formelle d'un schema de resolution de 1 'equation 
des ondes : I'erreur de methode 

Resume : Les schenias numcriqucs de resolution de I'cquation des ondes en dimension 1 sont 
notoirement convergents. Nous presentons una formalisation detaillee du schema le plus simple 
et nous prouvons formellement sa convergence dans le systeme d'aide a la preuve Coq. La difS- 
culte principale se situe dans la definition adequate des comportements asymptotiques et du fait 
que ces comportements ne sont pas completement decrits dans les preuves sur papier. A notre 
connaissance, c'est la premiere mecanisation complete d'une telle preuve en analyse numerique. 

Mots-cles : equations aux derivees partielles, equation des ondes acoustiques, schema numerique, 
preuve formelle en Coq 
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1 Introduction 

Ordinary differential equations (ODE) and partial differential equations (PDE) are ubiquitous in 
engineering and scientific computing. They show up in weather forecast, nuclear simulation, etc., 
and more generally in numerical simulation. Solutions to nontrivial problems are nonanalytic, 
hence approximated by numerical schemes over discrete grids. 

Numerical analysis is mainly interested in proving the convergence of these schemes, that 
is, the approximation quality increases as the size of the discretization steps decreases. The 
approximation quality is characterized by the error defined as the difference between the exact 
continuous solution and the approximated discrete solution; this error must tend toward zero in 
order for the numerical scheme to be useful. 

There is a wide literature on this topic, e.g. see [22l[23|, but no article goes into all the details. 
These "details" may have been skipped for readability, but they could also be mandatory details 
that were omitted due to an oversight. The purpose of a mechanically-checked proof is to uncover 
these issues and check whether they could jeopardize the correctness of the schemes. 

This work is a first step toward the development of formal tools for dealing with the convergence 
of numerical schemes. It would have been sensible to start with classical schemes for ODE, such 
as the Euler or Runge-Kutta methods. But we decided to directly validate the feasibility of 
our approach on the more complicated PDE. Moreover, this opens the door to a wide variety of 
applications, as they appear in many realistic problems from industry. 

We chose the domain of wave propagation because it represents one of the most common phys- 
ical phenomena one experiences in everyday life: directly through sight and hearing, but also via 
telecommunications, radar, medical imaging, etc. Industrial applications range from aeroacoustics 
to music acoustics (acoustic waves), from oil prospection to nondestructive testing (elastic waves), 
from optics to stealth technology (electromagnetic waves), and even include stabilization of ships 
and offshore platforms (surface gravity waves) . We restrained ourselves to the simplest example of 
wave propagation models, the acoustic wave equation in a one-dimensional space domain, for it is 
a prototype model for all other kinds of wave. In this case, the equation describes the propagation 
of pressure variations (or sound waves) in a fluid medium; it also models the behavior of a vibrat- 
ing string. For simplicity, we only consider homogeneous media, meaning that the propagation 
velocity is constant. Among the wide variety of numerical schemes for approximately solving the 
ID acoustic wave equation, we chose the simplest one: the second order centered finite difference 
scheme, also known as the "three-point scheme" . Again, for simplicity, we only consider regular 
grids with constant discretization steps for time and space. 

To our knowledge, this is the first time this kind of mathematical proof is machine-checked|^ 
Few works have been done on formalization and proofs on mathematical analysis inside proof 
assistants, and fewer on numerical analysis. Even real analysis developments are relatively new. 
The first developments on real numbers and real analysis are from the late 90 's [TUl [HI [TTl [T^ [T^ . 
Some intuitionist formalizations have been realized by a team at Nijmegen J13l l9]. Analysis results 
are available in provers such as ACL2, Coq, HOL Light, Isabelle, Mizar, or PVS. Regarding 
numerical analysis, we can cite [20] which deals, more precisely, with the formal proof of an 
automatic differentiation algorithm. About K" and the dot product, an extensive work has been 
done by Harrison [TS] . About the big O operator for asymptotic comparison, a decision procedure 
has been developed in [2] ; unfortunately, we needed a more powerful big O and those results were 
not applicable. 

Section [5] presents the PDE, the numerical scheme, and their mathematical properties. Sec- 
tion [3] describes the basic blocks of the formalization: dot product, big O, and Taylor expansions. 
Section m is devoted to the formal proof of the convergence of the numerical scheme. 



^ The Coq sources of the formal development are available from|http : //f ost . saclay ■ Inria . f r/Mave_method_error ■ php | 
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2 Wave Equation 

A partial differential equation modeling an evolutionary problem is an equation involving par- 
tial derivatives of an unknown function of several independent space and time variables. The 
uniqueness of the solution is obtained by imposing additional conditions, typically the value of the 
function and the value of some of its derivatives at the initial time. The right-hand sides of such 
initial conditions are also called Cauchy data, making the whole problem a Cauchy problem, or an 
initial-value problem. 

The mathematical theory is simpler when unbounded domains are considered |22j . When the 
space domain is bounded, the computation is simpler, but we have to take reflections at domain 
boundaries into account; this models a finite vibrating string fixed at both ends. Thanks to 
the nice property of finite velocity of propagation of the wave equation, we can build two Cauchy 
problems, one bounded and the other one unbounded, that coincide on the domain of the bounded 
one. Thus, we can benefit from the best of both worlds: the bounded problem makes computation 
simpler and the unbounded one avoids handling reflections. This section, as well as the steps taken 
at section 0] to conduct the proof of the convergence of the mmrerical scheme, is inspired by [3]. 

2.1 The continuous equation 

The chosen PDE models the propagation of waves along an ideal vibrating elastic string, see [BIT]. 
It is obtained from Newton's laws of motion [21] . 

The gravity is neglected, hence the string is supposed rectilinear when at rest. Let u{x,t) be 
the transverse displacement of the point of the string of abscissa x at time t from its equilibrium 
position. It is a (signed) scalar. Let c be the constant propagation velocity. It is a positive number 
that depends on the section and density of the string. Let s{x,t) be the external action on the 
point of abscissa x at time t; it is a source term, such that i = => s{x, t) = 0. Finally, let uo{x) 
and Ui{x) be the initial position and velocity of the point of abscissa x. We consider the Cauchy 
problem (i.e., with conditions at i = 0) 

1 f 3 u 
{L{c)u){x,t) = —{x,t) + A{c)u{x,t) = s{x,t), 

(Lim)(x,0) = —{x,<d)=ui{x), 
(Lqu){x,{)) — u{x,0) — Uq{x) 
where the differential operator A{c) is defined by 

(4) -^w"-^!.' 

We admit that under reasonable conditions on the Cauchy data uq and ui and on the source 
term s, there exists a unique solution to the Cauchy problem (IT])-(I3]) for each c > 0. This is a 
mathematical known fact (established for example from d'Alembcrt's formula ^), that is left 
unproved here. 

For such a solution u, it is natural to associate at each time t the positive definite quadratic 
quantity 



(1) 


Vi > 0, Vx e 


(2) 


Va;e 


(3) 


Vze 



(5) E{c){u){t)^^'\ 



du 
x^-(x,t) 



1 2 

+ -||a-i-^M(a;,t)||^(^) 



where (w,w) = J^v{x)w{x)dx, \\v\\ = {v,v) and ||w||^/^^ = {A{c)v,v). The first term is inter- 
preted as the kinetic energy, and the second term as the potential energy, making E the mechanical 
energy of the vibrating string. 
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Figure 1: Three-point scheme: u (x) depends on u^_^^ Uj, u^j^^ and u (•). 



This simple partial derivative equation happens to possess an analytical solution given by the 
so-called d'Alembert's formula [T7], obtained from the method of characteristics [TB], Vi > 0, 

1 -1 px + ct 

(6) u{x,t)^ -{u(i{x-'Ct) + U(i{x + ct)) + — I ui{y)dy+ 

2 2c 7^_ct 




a;4-c(i — (t) 



s(2/, (^)dy da. 



x — c{t — a) 



One can deduce from formula ([5]) the useful property of finite velocity of propagation. Assuming 
that we are only interested in the resolution of the Cauchy problem on a compact time interval of 
the form [0, tmax] with imax > 0, we suppose that uq, ui and s have a compact support. Then the 
property states that there exists Xmin and x^ax with Xmin < Xmax such that the support of the 
solution is a subset of J7 = [xmin, a^max] X [0, imax]- Furthermore, since the boundaries do not have 
time to be reached by the signal, the Cauchy problem set on H. by adding homogeneous Dirichlet 
boundary conditions {i.e. for all t £ [0,iinax], u{xuiin,t) = u(xniax,0 ~ 0)i admits the same 
solution. Hence, we will numerically solve the Cauchy problem on fl, but with the assumption 
that the spatial boundaries arc not reached. 

Note that the implementation of the compact spatial domain [a;min,a;max] will be abstracted 
by the notion of finite support (that is to say, being zero outside of an interval, see Section 14. 2p 
and will not appear explicitly otherwise. 

Note also that most properties of the continuous problem proved unnecessary in the formaliza- 
tion of the numerical scheme and the proof of its convergence. For instance, integration operators 
and d'Alembert's formula can be avoided as long as we suppose the existence and regularity of a 
solution to the PDF and that this solution has a finite support. 



2.2 The discrete equations 



Let (Ax, At) be a point in the interior of 17; define the discretization functions jax{x) = [ ^ ^°' 



main SI is approximated by the regular discrete grid defined by 



and kAtit) = [^j J ; then set jmax = JAxixma^x) and kn 



(). Now, the compact do- 



(7) Vke [0..fc,„ax],Vj€ [0..jn.ax], 

Let Vh be a discrete function over [0..ji, 



k dcf 



x:; = ixj,t^) 



feN dcf 



(a;min + jAx.kAt). 
X [0..fc,„ax]- For all k in [0..fcniax], we write v^ 



{^'j)o<j<jmB.^' then Vh — {{v^)o<k<k^^^) ■ A function v defined over fl is approximated at the points 

of the grid by the discrete function v^ defined on [0..j,nax] x [Cfemax] by v'- ~ v(x'^), except for u 

where we use the notation u^ ~ ''^{^j) to prevent notation clashes. 

Let uoh and uih be two discrete functions over [0..jniax]; let s/i be a discrete function over 
[0..jmax] X [0..fcniax]- Then, thc discrete function Uh over [0..jinax] x [0..fcmax] is said to be the 
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solution of the three-poinio finite difference scheme, as illustrated in Figure [TJ when the following 
set of equations holds: 



(8) Vfce [2..fc,nax],Vje [O.J-nax], 






(9) Vj e [O..J.„ax], {L,h(c) Uh), '^' ^^^^ + ^(A,(c) <), = ^i,„ 

(10) Vj e [0..j,„ax], {LohUh)] == w" = ItOj, 
(f f) Vfc e [0..fc,„ax[, W^l = wL, + i = 

where the matrix Ah{c), discrete analog of A(c), is defined, for any vector Vh = {{vj)o<j<k„^^y,), by 



(12) VjG[0..j„,ax], {Ah{c)vh)^ 



Ax 



,^2 



Note that defining Uh for artificial indexes j = — f and j = jmax + 1 is a trick to make the 
three-point spatial scheme valid for j — Q and j = jmax- 
A discrete analog of the energy is also defined bjo 



(13) E,,{c){uu) 2 



.2 

fe+i dcf 1 






M 



+ \{ulA'-') 



Ax 



2\^h,^h /A^{c) 



where {vh,Wh)^^ ^ J2^jTo ^j^j^^^ hhW^x == {vh,Vh)/^^, 

and {vh ,Wh)Ah (c) = ("^^ (c) Vh,'Wh)^^- 

Note that the three-point scheme is parametrized by the discrete Cauchy data mq/i and wi/i, and 
by the discrete source term Sh. Of course, when ug/i, ui^, and s;; are respectively approximations 
of Up, wi, /, then the discrete solution Uh is an approximation of the continuous solution u. 

2.3 Convergence 

Let ( and ^ be in ]0,1[ with ^ < 1 — ■C- The CFL(C, ^) condition (for Courant-Friedrichs-Lewy, 
see [5]) states that the discretization steps satisfy the relation 

(M) C.^.l-«^ 

Note that the lower bound (^ may seem surprising from a numerical analysis point of view; the 
formalization has however shown that it was mandatory (see Section |4?3|) . 

The convergence error Ch measures the distance between the continuous and discrete solutions, 
ft is defined by 

(15) Vfc e [0..fc,nax], Vj G [O..J„,ax], (^ =' t^ ~ u) . 

The truncation error Eh measures at which precision the continuous solution satisfies the numerical 
scheme. It is defined by 



(16) Vfc e [2..fc,nax], Vj G [O..J,nax], 4'' = ^^h{c) Uj^ " sf 

(17) Vj G [0..j,nax], 4 = {Lih{c)uh)j -Ul,j, 

(18) Vj G [0..j,nax], e~'^ ^ [LohUh)] -Ui3j. 



^In the sense "three spatial points", for the definition of matrix Ah(c). 

^By convention, the energy is defined between steps k and fc -|- 1, thus the notation k + ■^. 
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The numerical scheme is said to be convergent of order 2 if the convergence error tends toward 
zero at least as fast as Ax^ + Ai^ when both discretization steps tend toward 0. More precisely, 
the numerical scheme is said to be convergent of order (p,g) uniformly on the interval [0,iniax] if 
the convergence error satisfies (see Section [X^ for the definition of the big O notation that will be 
uniform with respect to space and time) 



(19) 



fcAt(t) 



o 



Ax 



[0J„ 



,(Aa;P + Ai«). 



The numerical scheme is said to be consistent with the continuous problem at order 2 if the 
truncation error tends toward zero at least as fast as Ax^ + At^ when the discretization steps 
tend toward 0. More precisely, the numerical scheme is said to be consistent with the continuous 
problem at order (p, q) uniformly on interval [0,tinax] if the truncation error satisfies 



(20) 



.feAt(t) 



0[o,t^_]iAxP + At^). 



The numerical scheme is said to be stable if the discrete solution of the associated homogeneous 
problem (i.e. without any source term, s(x,t) = 0) is bounded from above independently of the 
discretization steps. More precisely, the numerical scheme is said to be stable uniformly on interval 
[0, imax] if the discrete solution of the problem without any source term satisfies 



(21) 3a,Ci,C2>0, Vie [0,t„ 



, VAx, At > 

feAt(i) 



VAaJ+At2 < 



Ax 



< (Ci + C2t){\\uoh\\^^ + \\uoh\\A^{c) + hlhWAx)- 



The result to be formally proved at section |4] states that if the continuous solution u is regular 
enough on il and if the discretization steps satisfy the CFL(C, ^) condition, then the three-point 
scheme is convergent of order (2, 2) uniformly on interval [0,imax]- 

We do not admit (nor prove) the Lax equivalence theorem which stipulates that for a wide 
variety of problems and numerical schemes, consistency implies the equivalence between stability 
and convergence. Instead, we establish that consistency and stability implies convergence in the 
particular case of the one-dimensional acoustic wave equation. 

3 The Coq Formalization: Basic Blocks 

We decided to use the Coq proof assistant j4| , as Coq was already used to prove the floating-point 
error |S] of this case study. All our developments use the Coq real standard (classical) library. 
Numerical equations, numerical schemes, numerical approximations deal with classical statements, 
and are not in the scope of intuitionist theory. 



3.1 Dot product 

The function space Z — > R can be equipped with pointwise addition and multiplication by a 
scalar. The result is a vector space. In the following, we are only interested in functions with 
finite support, that is the subset 



i^='{/:Z 



3a, & e Z, Vi e Z, f{i) ^O^a<i<b}, 



which is also a vector space. Then it is possible to define a dot product on F, noted 
follows: 



as 



(22) 



(/,.9) = E/«3« 



dcf 



and the corresponding norm ||/|j = \/{f, /). The corresponding Coq formalization is not imme- 
diate, though. One could characterize F with a dependent type, but that would make operation 
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(., .) difficult to use (each time it is applied, proofs of finite support properties have to be passed 
as well). Instead, we define (.,.) on the full function space Z — > R using Hilbert's £-operator 
(provided in Coq standard library in module Epsilon), as follows: 

(23) (/, g) = e \ Xx3a b, (Vi, (/(*) ^ V .9(1) ^ 0) ^ a < z < 6) 

A2: = EL/«5(*) 

Said otherwise, wc give (/, g) a definition as a finite sum whenever / and g both have finite support 
and we let (/, g) undefined otherwise. 

To ease the manipulation of functions with finite support, we introduce the following predicate 
characterizing such functions 

FS{f) =*' 3a6,Vi, f(i) ^^ ^ a < i < b 

and we prove several lemmas about it, such as 

Vfg,FS{f)^FS{g)^FS{f + g) 
yfc, FSif) => FS{c ■ f) 
yfk,FS{f)^FS{i^f{i + k)) 

We also provide a Coq tactic to automatically discharge most goals about FS{.). Finally, we can 
establish lemmas about the dot product, provided functions have finite support. Here are some of 
these lemmas: 

V/ g c, FS{.f) => FS{g) => (c • /, g)=c- (/, g) 

V/i h g, FS{h) => FS{f2) ^ FS{g) => {h + .f2, ,9) = (./i, ,9) + {h, g) 

\Jfg,FS{f)^FS{g)^\{f,g)\<\\f\\-\\g\\ 

yfg,FS{f) ^ FS{g) ^ ||/ + .9|| < ||/|| + ||g|| 

These lemmas are proved by reduction to finite sums, thanks to Formula ([25)1 . Note that the value 
of {f,g)^.^ defined in Section [2?2] is equal to Ax ■ {f,g)- 

3.2 Big O notation 

For two functions / and g over M", one usually writes f{x) = Op||^o(9(^)) for 
3a, 00, VfeM", ||.f|| <a=> |/(:e)| <C-|.9(f)|. 

Unfortunately, this definition is not sufficient for our formalism. Indeed, while /(x. Ax) will 
be defined over M^ x M^, ^(Ax) will be defined over M^ only. So it begs the question: what to do 
about X? 

Our first approach was to use 

Vx, /(x,Ax) = 0||Ax|Ho(g(Ax)) 

that is to say 

Vx,3a,C>0, VAxeM^ || Ax|| < a =^ |/(x, Ax)| < C • |g(Ax)| 

which means that a and C are functions of x. So we would need to take the minimum of all the 
possible values of a, and the maximum for C. Potentially, they may be and +00 respectively, 
making them useless. 

In order to solve this issue, we had to define a notion of big O uniform with respect to the 
additional variable x: 

3a, 00, Vx,Ax, ||Ax||<a^|/(x,Ax)|<C-|(7(Ax)|. 
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Variables x and Ax are restricted to subsets S and P of M^. For instance, the big O that 
appears in Equation (J19p uses 

S = Kx[0,t„,ax], 

P = i Ax = (Ax, Ai) I < Ax A < At A C < '-r-^ < 1 - ^ 

As often, the formal specification has allowed us to detect some flaws in usual mathematical 
pen-and-papcr proofs, such as an erroneous switching of the universal and existential quantifiers 
hidden in the big O definition. 

3.3 Taylor expansion 

The formalization assumes that "sufficiently regular" functions can be uniformly approximated by 
multivariate Taylor series. More precisely, the development starts by assuming that there exists 
two operators partial_derive Jirstvar and _secondvar. Given a real-valued function / defined 
on the 2D plane and a point of it, they respectively return the functions -^ and -J^ for this point, 
if they exist. 

Again, these operators are similar to the use of Hilbert's e operator. For documentation 
purpose, one could add two axioms stating that the returned function computes the derivatives 
for derivable functions; they are not needed for the later development though. Indeed, none of our 
proofs depend on the actual properties of derivatives; they only care about the fact that differential 
operators appear in both the regularity definition below and the wave equation. 

The two primitive operators -4- and ^ are encompassed in a generalized differential operator 
dx"^dt" • This allows us to define the 2D Taylor expansion of order n of a function /: 

DL,.(/,x) ^^' (A,io « t;^ (E (:) ■ s;5?-^(-) ■ ^-- ■ ^"-■) ■ 

p=Q ^ \m=0 ^ / / 

A function / is then said to be sufficiently regular of order n if 

(24) Vm<n, DL™_i(/, x)(Ax) - /(x + Ax) = O (|| Ax|r ) . 

4 The Coq Formalization: Convergence 

4.1 Wave equation 

As explained in Section [51 a solution of the wave equation with given uq, ui and s verifies Equa- 
tions (H}-®. Its discrete approximation verifies Equations ©-([TU]). Both are directly translated 
in Coq using the definitions of Section [31 Concerning the discretization, we choose that the space 
index is in Z (to be coherent with the dot product definition of Section 13. ip while the time index 
is in N. 

Our goal is to prove the uniform convergence of the scheme with order (2,2) on the interval 

U , trnax ■ 



fcAt(t) 



^tefcwxl (Ax2+At2) 



t e [0, tn,ax] 
(Ax, At) -> 

< Ax A < At A 
C<cfi<l-« 

4.2 Finite support 

The proofs concerning the convergence of the scheme rely on the dot product. As explained in 
Section [3?T1 the dot product requires the functions to have a finite support in order to apply any 
lemma. We therefore proved the finiteness of the support of many functions. We assume that the 
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inputs uq, ui, and s of the wave equation have a finite support. More precisely, we assume that 
there exists xi and X2 such that uo{x) = ui{x) = for all x out of [xi, X2] and s{x, t) = for all 
X out of [xi — c ■ t,X2 + c ■ t] where c is the velocity of propagation of waves in Equation ([1]). 

Figure [2] describes the nullity, that is to say the finite support, of the various functions. We 
needed to prove the finiteness of their support: 



I Mo and ui may be nonzero. 
^ s and thus u may be nonzero. 
V Uh may be nonzero. 

slope: H: • [c • |i] "' (equals || under CFL 
X conditions) 




Xi X2 



Figure 2: Finite supports. The support of the Cauchy data uq and ui is included in the support of 
the continuous source term s, and of the continuous solution u. Which is in turn also included in 
the support of the discrete solution m^, provided that the CFL condition holds. For a finite tmax, 
all these supports are finite. 

• Mo and Ml by hypothesis and therefore uqj and mi.j. 

• s (for any value t) by hypothesis and therefore s^ is zero outside of a cone of slope c~^. 

• the scheme itself has a finite support: due to the definition of m^ and the nullity of mq.j and 
Ml J- and s^, we can prove that m^' is zero outside of a cone of slope -^ • [c • ■^~\ . Under 
CFL(C,'^) conditions, this slope will be ■^. 

• the truncation and convergence errors also have finite support with the previous slope. 

We need here an axiom about the nullity of the continuous solution. We assume that the 
continuous solution m(x, t) is zero for x out of [xi — c-i,X2+c-t] (same ass). This is mathematically 
correct, since it derives from d'Alembert's formula ([6]). But its proof is out of the scope of the 
current formalization and wc therefore preferred to simply add the nullity axiom. 

4.3 Consistency 

We first prove that the truncation error is of order Ax^ + At^. The idea is to show that, for Ax 
small enough, the values of the scheme L^ are near the corresponding values of L. This is done 
using the properties of Taylor expansions. This involves long and complex expressions but the 
proof is straightforward. 

We first prove that the truncation error in one point (j, k) is a 0{Ax^ + Ai^). This is proved 
for fc = and fc = 1 by taking advantage of the initializations and Taylor expansions. For bigger fc, 
the truncation error reduces to the sum of two Taylor expansions of degree 3 in time (this means 
771 = 4 in Formula (j24p ) and two Taylor expansions of degree 3 in space that partially cancel 
(divided by something proportional to || Ax|p). Here, we take advantage of the generality of big 
O as we consider the sum of a Taylor expansion on Ax and of a Taylor expansion on —Ax. If we 
had required < Ax (as a space grid step), we could not have done this proof. 

The most interesting part is to go from pointwise consistency to uniform consistency. We want 
to prove that the norm of the truncation error (in the sense of the infinite dot product {■,■) ^^) is 
also 0{Ax^ + At^). We therefore need to bound the number of nonzero values of the truncation 
error. As explained in Section 14.21 the truncation error values at time k ■ At may be nonzero 
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between xi'k = LfeJ ~ T'' ' Sill ^ ^^'^ ^2^ = [^] + \c ■ -^J /c. This gives a number of terms N 
roughly bounded by (aU details are handled in the formal proof): 



jY < X2fc " Xifc ^ X2 



< 



Xi 



Ax 

X2 -Xi 
Ax2 



Aa;2 



Z ■ /CiT 



At 



2- 



Ax 



Ai 



Ax 



Ax 



As the norm is a Ax-norm, this reduces to bounding with a constant value the value N ■ Ax^ 



which is smaller than X2 — Xi + 2 • imax ■ c + 2 ■ t„ 



Ax 
At ■ 



To bound this with a constant value, 



we require c-^ to have a constant lower bound C (it already had an upper bound 1 — ^). Then 



iV • Ax^ < X2 - XI + 2 • i„ 



2 • c • t 



max ■ 7 which is constant. 



Mathematically, this requirement comes as a surprise. The following scenario explains it. If 
goes to zero, then At goes to zero much faster than Ax. It corresponds to Figure [31 The 
number of nonzero terms (for Uh and thus for the truncation error) goes to infinity as -^ goes to 
zero. 



-A* 

'Ai 




?Ai 



Figure 3: For a given time to, the number of nonzero values increases when the slope -^ goes to 
zero. From left to right. At is divided by 2 whereas Ax remains the same. We can see that the 
number of nonzero terms is almost doubled (from 9 to 17). 



4.4 Stability 

To prove stability, we use the discrete energy defined in Equation ([TS]) . From the properties of the 
scheme, we calculate the evolution of the energy. At each step, it increases by a known value. In 
particular, if s is zero, the discrete energy (as the continuous energy) is constant: 

Vfc > 0, E,ic)iu^f+i - E,ic)iu,f-i = i («|;+i - ut\ sf;)^, . 



From this, we give an underestimation of the energy: 



V/c, 




,k+l 



At 



< Eh{c){uh) 



k+i 



Therefore we have the nonnegativity of the energy under CFL(C,^) conditions. For convergence, 
the key result is the overestimation of the energy: 

V2 



^ Eh{c){uhf+-- < ^ Eh{c){uh)-^ + 



2^2^^ ^t■^¥^s,i^,mA. 



for all time i, with k = [^J — 1. 

This completes the stability proof. In the inequality above, the energy is bounded for u/j. but 
the bound is actually valid for all the solutions of the discrete scheme, for any initial conditions 
and source term. 

Note that the formal proof of stability closely follows the mathematical pen-and-paper proof 
and no additional hypotheses were found to be necessary. 
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4.5 Convergence 

Wc prove that the convergence error is the solution of a scheme and therefore the results of 
Section 23] apply to it. More precisely, for all Ax. the convergence error is solution of a discrete 
scheme with inputs 

uq.j = 0, Ml J = ■^, and s)' = £^^+\ 

where the errors refer to the errors of the initial scheme of the wave equation with grid steps Ax. 
(Actual Coq notations depend on many more variables.) 

We have proved many lemmas about the initializations of our scheme and of the convergence 
error. The idea is to prove that the initializations of the scheme are precise enough to guarantee 
that the initial convergence errors (at step and 1) are accurate enough. 

We also bounded the energy of the convergence error. Using results of Section 231 the proof 
reduces to bounding the sum of the source terms, here the truncation errors. Using results of 
Section l473l we prove this sum to be 0{Ax^ + At^). A few more steps conclude the proof. 

Once more, the formal proof follows the pen-and-paper proof and progresses smoothly under 
the required hypothesis, including all the conditions on ^ of Equation ([T^. 

5 Conclusion and perspectives 

One of the goals of this work is to favor the use of formal methods in numerical analysis. It may 
seem to be just wishful thinking, but it is actually seen as needed by some applied mathematicians. 
An early case led to the certification of the O^yssee tool [5D|. This tool performs automatic 
differentiation, which is one of the basic blocks for gradient-hased algorithms. Our work tackles 
the converse problem: instead of considering derivation-based algorithms, we have formalized and 
proved part of the mathematical background behind integration-based algorithms. 

This work shows there may be a synergy between applied mathematicians and logicians. Both 
domains are required here: applied mathematics for an initial proof that could be enriched upon 
request and formal methods for machine-checking it. This may be the reason why such proofs are 
scarce as this kind of collaboration is uncommon. 

Proof assistants seem to mainly deal with algebra, but we have demonstrated that formalizing 
numerical analysis is possible too. We can confirm that pen-and-paper proofs are sometimes 
sketchy: they may be fuzzy about the needed hypotheses, especially when switching quantifiers. 
We have also learned that filling the gaps may cause us to go back to the drawing board and to 
change the basic blocks of our formalization to make them more generic (a big O that needs to be 
uniform and also generic with respect to a property P). 

The formal bound on the error method, while of mathematical interest, is not sufficient to 
guarantee the correction of numerical applications implementing the three-point scheme. Indeed, 
such applications usually perform approximated computations, e.g., floating-point computations, 
for efficiency and simplicity reasons. As a consequence, the proof of the method error has to 
be combined with a proof on the rounding error, in order to get a full-fledged correction proof. 
Fortunately, the proof on the rounding error has already been achieved [5] • We are therefore close 
to having a formal proof of both the numerical scheme and its floating-point implementation. 

An advantage of Coq with respect to most other proof assistants is the ability to extract pro- 
grams from proofs [18j . For this work, it does not make much sense to extract the algorithm from 
the proofs: not only is the algorithm already well-known, but its floating-point implementation 
was also certifled [6]. So, an extraction of the algorithm would not bring much. However, ex- 
traction gives access to the constant C hidden behind the big O notation. Indeed, the proof of 
the floating-point algorithm relies on the discrete solution being good enough, so that the com- 
puted result does not diverge. Precisely, the convergence error has to be smaller than 1, and an 
extracted computation would be able to ensure this property. Furthermore, having access to this 
constant can be useful to the applied mathematicians for the a posteriori estimations needed for 
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adaptive mesh refinements. Extraction also gives access to the a constant. That way, we could 
check that the constant Ax chosen in the C program described in [5] verifies this requirement. 
Note that performing an extraction requires to modify the definition of the big O so that it lives in 
Set instead of Prop. But this formalization change happens to be straightforward and Coq then 
succeeds in extracting mathematical formulas for constants a and C. Only basic operators {e.g. 
+, v^, min) and constants {e.g. <max, £,^ Xii Taylor constants) appear in them, so they should be 
usable in practice. 

The formal development is about 4500-linc long. Its dependency graph is detailed in Figure SI 
About half of the development is a reusable library described in Section |2] and the other half is the 
proof of convergence of the numerical scheme described in Section H) This may seem a long proof 
for a single scheme for a single PDE. To put it into perspective, usual pcn-and-paper proofs are 
10-page long and an in-depth proof can be 60-page long. (We wrote one to ensure that we were 
not getting sidetracked.) So, at least from a length point of view, the formal proof is comparable 
to a detailed pen-and-paper proof. 
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Figure 4: Dependency graph of the Coq development. On the left are the files from the convergence 
proof. The other files correspond to the reusable library. 

In the end, the whole development contains only two axioms: the e operator for the infinite dot 
product (see Section 13. ip and the finite support of the continuous solution of the wave equation 
(see Section [4.2^ . So, except for this last axiom which is related to the chosen PDE, the full 
numerical analysis proof of convergence is machine- checked and all required hypotheses are made 
clear. There is no loss of confidence due to this axiom, since the kind of proof and the results it 
is based upon are completely different from the ones presented here. Indeed, this axiom is about 
continuous solutions and hence much less error-prone. 

For this exploratory work, we only considered the simple three-point scheme for the one- 
dimensional wave equation. Further works involve generalizing our approach to other schemes 
and other PDEs. We are confident that it would scale to higher-dimension and higher-order equa- 
tions solved by discrete numerical schemes. However, the proofs of Section U are entangled with 
particulars of the presented problem, and would therefore have to be redone for other problems. 
So a more fruitful approach would be to prove once and for all the Lax equivalence theorem that 
states that consistency implies the equivalence between convergence and stability. This would 
considerably reduce the amount of work needed for tackling other schemes and equations. 

This work also showed us that summations and finite support functions play a much more 
important role in the development than we first expected. We are therefore considering moving 
to the SSReflect interface and libraries for Coq [S], so as to simplify the manipulations of these 
objects in our forthcoming works. 
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